/**
 * Copyright 2009-2010 - YangJiandong(chunquedong)
 * 
 * This file is part of ChunMap project
 * Licensed under the GNU LESSER GENERAL PUBLIC LICENSE(Version >=3)
 * 你可以自由复制、传播本项目的下载包文件，但必须保持其完整性。
 * 我们不用对使用中的风险及由此造成的损失承担任何责任。
 * 详细情况请见《ChunMap许可协议》。

 * 想了解更多有关ChunMap的信息，请访问http://code.google.com/p/chunmap/
 */
package chunmap.util.math;

import java.util.Arrays;

import chunmap.util.ChunmapException;

/**
 * 矩阵类
 * @author chunquedong
 * 
 */
public class Matrix implements Cloneable {
	private double[][] a;
	// m行n列
	private final int m, n;
	private String name;

	public Matrix(int am, int an) {
		m = am;
		n = an;
		a = new double[m][n];
		name = "matrix" + m + "*" + n;
	}

	public Matrix(int am, int an, String aName) {
		m = am;
		n = an;
		a = new double[m][n];
		name = aName;
	}

	public int getM() {
		return m;
	}

	public int getN() {
		return n;
	}

	public void set(int r, int c, double value) {
		a[r][c] = value;
	}

	public double get(int r, int c) {
		return a[r][c];
	}

	@Override
	public String toString() {
		String s;
		s = name + ":\n";

		int m = getM();
		int n = getN();

		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				s += String.format("%f", get(i, j)) + "\t";
			}
			s += "\n";
		}
		return s;
	}

	public String getName() {
		return name;
	}

	public void setName(String name) {
		this.name = name;
	}

	@Override
	public Object clone() {
		Object obj = null;
		double[][] d = new double[m][n];
		try {
			obj = super.clone();
		} catch (CloneNotSupportedException e) {
			e.printStackTrace();
		}
		Matrix b = (Matrix) obj;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				d[i][j] = this.get(i, j);
			}
		}
		b.a = d;
		b.setName(name + "copy");
		return b;
	}

	@Override
	public int hashCode() {
		final int prime = 31;
		int result = 1;
		result = prime * result + Arrays.hashCode(a);
		result = prime * result + m;
		result = prime * result + n;
		return result;
	}

	@Override
	public boolean equals(Object obj) {
		if (this == obj)
			return true;
		if (obj == null)
			return false;
		if (getClass() != obj.getClass())
			return false;
		Matrix other = (Matrix) obj;
		if (!equals(this.a, other.a))
			return false;
		if (m != other.m)
			return false;
		if (n != other.n)
			return false;
		return true;
	}

	private boolean equals(double[][] a, double[][] b) {
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < m; j++) {
				if (!MyDouble.approximateEquals(a[i][j], b[i][j]))
					return false;
			}
		}
		return true;
	}

	public static Matrix unitMatrix(int n) {
		Matrix b = new Matrix(n, n, n + "unit");
		for (int i = 0; i < n; i++) {
			b.set(i, i, 1);
		}
		return b;
	}

	// ------------------------------------------------------------
	// 以下代码是我以前在c#上写的
	// ------------------------------------------------------------
	/**
	 * 矩阵转置
	 */
	public Matrix matrixTrans() {
		Matrix c = new Matrix(n, m);

		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				c.set(i, j, get(j, i));

		return c;
	}

	/**
	 * 矩阵加法
	 */
	public Matrix matrixAdd(Matrix b) {
		int m2 = b.getM();
		int n2 = b.getN();

		if ((m != m2) || (n != n2)) {
			throw new IllegalArgumentException("数组维数不匹配");
		}

		Matrix c = new Matrix(m, n);

		int i, j;
		for (i = 0; i < m; i++)
			for (j = 0; j < n; j++)
				c.set(i, j, get(i, j) + b.get(i, j));

		return c;
	}

	/**
	 * 矩阵数乘
	 */
	public Matrix matrixMulti(double k) {
		Matrix c = new Matrix(m, n);

		int i, j;
		for (i = 0; i < m; i++)
			for (j = 0; j < n; j++)
				c.set(i, j, get(i, j) * k);

		return c;
	}

	/**
	 * 矩阵乘法
	 */
	public Matrix matrixMulti(Matrix b) {
		int m2 = b.getM();
		int n2 = b.getN();

		if (n != m2) {
			throw new IllegalArgumentException("数组维数不匹配");
		}

		Matrix c = new Matrix(m, n2);

		int i, j, k;
		for (i = 0; i < m; i++) {
			for (j = 0; j < n2; j++) {
				double value = 0;
				for (k = 0; k < n; k++) {
					value += get(i, k) * b.get(k, j);
				}
				c.set(i, j, value);
			}
		}

		return c;
	}

	/**
	 * 逆矩阵,算法：高斯法
	 * 
	 * @return
	 * @throws MatrixException 
	 */
	public Matrix matrixInv() throws MatrixException {
		if (m != n) {
			throw new IllegalArgumentException("数组维数不匹配");
		}
		Matrix b = new Matrix(m, n);
		Matrix a = (Matrix) this.clone();

		int i, j, row, k;
		double max, temp;

		// 单位矩阵
		for (i = 0; i < n; i++) {
			b.set(i, i, 1);
		}
		for (k = 0; k < n; k++) {
			max = 0;
			row = k;
			// 找最大元，其所在行为row
			for (i = k; i < n; i++) {
				temp = Math.abs(a.get(i, k));
				if (max < temp) {
					max = temp;
					row = i;
				}

			}
			if (max == 0) {
				throw new MatrixException("没有逆矩阵");
			}
			// 交换k与row行
			if (row != k) {
				for (j = 0; j < n; j++) {
					temp = a.get(row, j);
					a.set(row, j, a.get(k, j));
					a.set(k, j, temp);

					temp = b.get(row, j);
					b.set(row, j, b.get(k, j));
					b.set(k, j, temp);
				}

			}

			// 首元化为1
			for (j = k + 1; j < n; j++)
				a.set(k, j, (a.get(k, j) / a.get(k, k)));
			for (j = 0; j < n; j++)
				b.set(k, j, (b.get(k, j) / a.get(k, k)));

			a.set(k, k, 1);

			// k列化为0
			// 对a
			for (j = k + 1; j < n; j++) {
				for (i = 0; i < k; i++)
					a.set(i, j, a.get(i, j) - (a.get(i, k) * a.get(k, j)));
				for (i = k + 1; i < n; i++)
					a.set(i, j, a.get(i, j) - (a.get(i, k) * a.get(k, j)));
			}
			// 对b
			for (j = 0; j < n; j++) {
				for (i = 0; i < k; i++)
					b.set(i, j, b.get(i, j) - (a.get(i, k) * b.get(k, j)));
				for (i = k + 1; i < n; i++)
					b.set(i, j, b.get(i, j) - (a.get(i, k) * b.get(k, j)));
			}
			for (i = 0; i < n; i++)
				a.set(i, k, 0);
			a.set(k, k, 1);
		}

		return b;
	}

	/**
	 * 矩阵求逆(伴随矩阵法)
	 */
	public Matrix matrixInvByCom() {
		double d = matrixDet();
		if (d == 0) {
			throw new IllegalArgumentException("没有逆矩阵");
		}
		Matrix Ax = matrixCom();
		Matrix An = Ax.matrixMulti(1.0 / d);
		return An;
	}

	/**
	 * 对应行列的代数余子式矩阵
	 */
	public Matrix matrixSpa(int ai, int aj) {
		if (m != n) {
			throw new IllegalArgumentException("数组维数不匹配");
		}
		int n2 = n - 1;
		Matrix b = new Matrix(n2, n2);

		// 左上
		for (int i = 0; i < ai; i++)
			for (int j = 0; j < aj; j++) {
				b.set(i, j, get(i, j));
			}
		// 右下
		for (int i = ai; i < n2; i++)
			for (int j = aj; j < n2; j++) {
				b.set(i, j, get(i + 1, j + 1));
			}
		// 右上
		for (int i = 0; i < ai; i++)
			for (int j = aj; j < n2; j++) {
				b.set(i, j, get(i, j + 1));
			}
		// 左下
		for (int i = ai; i < n2; i++)
			for (int j = 0; j < aj; j++) {
				b.set(i, j, get(i + 1, j));
			}
		// 符号位
		if ((ai + aj) % 2 != 0) {
			for (int i = 0; i < n2; i++)
				b.set(i, 0, -b.get(i, 0));

		}
		return b;
	}

	/**
	 * 矩阵的行列式
	 */
	public double matrixDet() {
		if (m != n) {
			throw new IllegalArgumentException("数组维数不匹配");
		}
		if (n == 1)
			return get(0, 0);

		double d = 0;
		for (int i = 0; i < n; i++) {
			d += get(1, i) * matrixSpa(1, i).matrixDet();
		}
		return d;
	}

	/**
	 * 矩阵的伴随矩阵
	 */
	public Matrix matrixCom() {
		Matrix c = new Matrix(m, n);

		for (int i = 0; i < m; i++)
			for (int j = 0; j < n; j++)
				c.set(i, j, matrixSpa(j, i).matrixDet());

		return c;
	}
	
	public class MatrixException extends Exception implements ChunmapException{
		private static final long serialVersionUID = -1578205432238297337L;

		public MatrixException(String arg0) {
			super(arg0);
		}
	}
}
